This Jupyter notebook contains code developed by N. Swanson-Hysell and Y. Park that conducts analysis and develops visualizations associated with the following Report published in Science in March 2019:
Arc-continent collisions in the tropics set Earth’s climate state
Macdonald, F.M., Swanson-Hysell, N.L., Park, Y., Lisiecki, L., and Jagoutz, O.
This code is written in Python 2.7 using the following python modules that can be installed via pip or conda:
The full computational environment associated with these modules can be found within the suture_analysis.yml file in this repository. In addition, the code imports and uses custom python functions that are included within this repository as recon_tools.py. The implementation of this analysis relies heavily on functions within pygplates http://www.gplates.org/docs/pygplates/index.html. The pygplates code is compatible with Python 2, which is why this code is in Python 2 rather than 3.
import pygplates
import pandas as pd
import numpy as np
import random
import matplotlib.pyplot as plt
from matplotlib.legend import Legend
import matplotlib.patches as patches
import cartopy
import cartopy.crs as ccrs
from shapely.geometry.polygon import Polygon
from recon_tools import *
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
The following code imports the following files:
The different suture files are generated by opening the Suture_Line.shp file in GPlates with different Suture_Lines.shp.gplates.xml files and then saving as .gpml files.
sutures_CEED_max = pygplates.FeatureCollection.read('data/sutures/CEED_Ex_max/Suture_Lines.gpml')
sutures_CEED_max_min = pygplates.FeatureCollection.read('data/sutures/CEED_Ex_max_min/Suture_Lines.gpml')
sutures_M_max = pygplates.FeatureCollection.read('data/sutures/CEED_Ex_max/Suture_Lines.gpml')
sutures_M_max_min = pygplates.FeatureCollection.read('data/sutures/CEED_Ex_max_min/Suture_Lines.gpml')
CEED_blocks = pygplates.FeatureCollection.read('paleogeo_models/CEED6+Kazakh/CEED6+D18_Kazakh.shp')
CEED_land = pygplates.FeatureCollection.read('paleogeo_models/CEED6/CEED6_LAND.shp')
Matthews_blocks = pygplates.FeatureCollection.read('paleogeo_models/Matthews_etal_2016_Global_Plate_Model_GPC/ContPolys/PresentDay_ContinentalPolygons_Matthews++.shp')
TC17_SHM17_D18_model = 'paleogeo_models/TC2017_SHM2017_D2018.rot'
M2016_model = 'paleogeo_models/Matthews_etal_2016_Global_Plate_Model_GPC/Global_EB_250-0Ma_GK07_Matthews++.rot'
There are additional attributes associated with the ophiolite-bearing suture shapefile that cannot be mapped into GPlates fields. To be able to deal with and interrogate these, we build up a pandas Dataframe and a dictionary from a .csv filed exported from QGIS. The attributes for a suture can be accessed using the suture name.
suture_attributes_df = pd.read_csv('data/sutures/Suture_Lines_attributes.csv',index_col='NAME')
suture_attributes = suture_attributes_df.to_dict('index')
suture_attributes['Kunlun']
# input parameters
t = 0 #reconstruction time
anchor = 1 #anchor plate ID (1 : spin axis)
# reconstruct the continents to the present day
# (which seems a bit silly, but exposes features associated with the reconstructed object)
CEED_blocks_present = []
pygplates.reconstruct(CEED_blocks, TC17_SHM17_D18_model, CEED_blocks_present, t, anchor)
CEED_land_present = []
pygplates.reconstruct(CEED_land, TC17_SHM17_D18_model, CEED_land_present, t, anchor)
# carry out the present-day reconstruction for the sutures
sutures_present = []
pygplates.reconstruct(sutures_CEED_max, TC17_SHM17_D18_model, sutures_present, t, anchor)
# carry out the present-day reconstruction for the sutures
sutures_present_active = []
pygplates.reconstruct(sutures_CEED_max_min, TC17_SHM17_D18_model, sutures_present_active, t, anchor)
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
plot_reconstructed_feature(ax,CEED_blocks_present)
plot_reconstructed_feature(ax,sutures_present,color='green',linewidth=2)
plot_reconstructed_feature(ax,sutures_present_active,color='orange',linewidth=2.5)
plt.title('ophiolite-bearing sutures (orange are active at present-day)')
plt.show()
sutures_Cenozoic_collection = pygplates.FeatureCollection.read('data/sutures/Age_grouped/Suture_Lines_Cenozoic.shp')
sutures_Cenozoic = []
pygplates.reconstruct(sutures_Cenozoic_collection, TC17_SHM17_D18_model, sutures_Cenozoic, t, anchor)
sutures_Mesozoic_collection = pygplates.FeatureCollection.read('data/sutures/Age_grouped/Suture_Lines_Mesozoic.shp')
sutures_Mesozoic = []
pygplates.reconstruct(sutures_Mesozoic_collection, TC17_SHM17_D18_model, sutures_Mesozoic, t, anchor)
sutures_Paleozoic_late_collection = pygplates.FeatureCollection.read('data/sutures/Age_grouped/Suture_Lines_Paleozoic_late.shp')
sutures_Paleozoic_late = []
pygplates.reconstruct(sutures_Paleozoic_late_collection, TC17_SHM17_D18_model, sutures_Paleozoic_late, t, anchor)
sutures_Paleozoic_early_collection = pygplates.FeatureCollection.read('data/sutures/Age_grouped/Suture_Lines_Paleozoic_early.shp')
sutures_Paleozoic_early = []
pygplates.reconstruct(sutures_Paleozoic_early_collection, TC17_SHM17_D18_model, sutures_Paleozoic_early, t, anchor)
Here is a 20 palette Viridis color scheme. Let's grab equally spaced colors for suture coloring.
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
plot_reconstructed_feature(ax,CEED_blocks_present)
plot_reconstructed_feature(ax,sutures_Cenozoic,color='#440154',linewidth=2)
plot_reconstructed_feature(ax,sutures_Mesozoic,color='#33638D',linewidth=2)
plot_reconstructed_feature(ax,sutures_Paleozoic_late,color='#3CBB75',linewidth=2)
plot_reconstructed_feature(ax,sutures_Paleozoic_early,color='#FDE725',linewidth=2)
Cenozoic_patch = patches.Patch(color='#440154',linestyle='-', label='Cenozoic (0 to 66 Ma)')
Mesozoic_patch = patches.Patch(color='#33638D', label='Mesozoic (66 to 252 Ma)')
Paleozoic_late_patch = patches.Patch(color='#3CBB75', label='Paleozoic (250 to 440 Ma)')
Paleozoic_early_patch = patches.Patch(color='#FDE725', label='Paleozoic (440 to 541 Ma)')
plt.legend(ncol=2,bbox_to_anchor=(0.5, -.02),handles=[Cenozoic_patch,Mesozoic_patch,Paleozoic_late_patch,Paleozoic_early_patch],loc='lower center')
plt.show()
The calculated value below for the sum of the arc lengths is quite close to what one gets when calculating line lengths in QGIS for this shapefile.
suture_lengths_all = []
suture_lengths_active = []
for j in range(len(sutures_present)):
arc_length = sutures_present[j].get_reconstructed_geometry().get_arc_length()
arc_length_km = arc_length*pygplates.Earth.mean_radius_in_kms
suture_lengths_all.append(arc_length_km)
for j in range(len(sutures_present_active)):
arc_length = sutures_present_active[j].get_reconstructed_geometry().get_arc_length()
arc_length_km = arc_length*pygplates.Earth.mean_radius_in_kms
suture_lengths_active.append(arc_length_km)
print('total suture length calculate using pyGPlates')
print(np.sum(suture_lengths_all))
print('active suture length')
print(np.sum(suture_lengths_active))
print('total suture length calculated in QGIS')
print(suture_attributes_df.LENGTH_KM.sum())
The function get_lengths_in_bands within recon_tools.py splits the polylines based on the latitude minimums and maximums provided. Let's do the split on the active sutures and plot the resulting lines to visualize this split.
# create the latitude bands
band_width = 10
lat_mins = np.arange(-90 , 90 , band_width)
lat_maxs = np.arange(-90+band_width, 90+band_width, band_width)
lat_mids = lat_mins + (lat_maxs-lat_mins)/2
# perform the calculations for sutures
active0_lengths, active0_polylines, active0_attributes = get_lengths_in_bands(sutures_present_active, lat_mins, lat_maxs)
np.sum(active0_lengths)
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,10),linestyle='--')
plot_reconstructed_feature(ax,sutures_present_active,color='orange',linewidth=3)
for i in range(len(active0_polylines)):
for j in range(len(active0_polylines[i])):
random_color = np.random.rand(3)
# pull out lat/lon vertices
lat_lon_array = active0_polylines[i][j].to_lat_lon_array()
lats = lat_lon_array[:,0]
lons = lat_lon_array[:,1]
ax.plot(lons,lats, transform=ccrs.Geodetic(), color=random_color, linewidth=1.5)
plt.show()
With the lengths split into latitude bands, the latitudinal distribution of modern active sutures can be plotted as a histogram.
plt.bar(lat_mids,active0_lengths,width=8,align='center',color='darkgreen')
plt.xlim(-90,90)
plt.ylabel('suture length in latitude band (km)')
plt.xlabel('latitude')
plt.title('all active ophiolite-bearing sutures')
plt.show()
These times correspond to the peaks in the total suture length and are included in Figure 1 of the main text.
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
ax.add_patch(patches.Rectangle(xy=[-180, -15], width=360, height=30,
facecolor='#52BE80',
alpha=0.2,
transform=ccrs.PlateCarree()))
plot_reconstructed_feature(ax,CEED_land_present,color='#D7DBDD',linewidth=0.5)
plot_reconstructed_feature(ax,CEED_blocks_present,color='#909497')
plot_reconstructed_feature(ax,sutures_Cenozoic,color='#440154',linewidth=2)
plot_reconstructed_feature(ax,sutures_Mesozoic,color='#33638D',linewidth=2)
plot_reconstructed_feature(ax,sutures_Paleozoic_late,color='#3CBB75',linewidth=2)
plot_reconstructed_feature(ax,sutures_Paleozoic_early,color='#FDE725',linewidth=2)
Cenozoic_patch = patches.Patch(color='#440154',linestyle='-', label='Cenozoic (0 to 66 Ma)')
Mesozoic_patch = patches.Patch(color='#33638D', label='Mesozoic (66 to 252 Ma)')
Paleozoic_late_patch = patches.Patch(color='#3CBB75', label='Paleozoic (250 to 440 Ma)')
Paleozoic_early_patch = patches.Patch(color='#FDE725', label='Paleozoic (440 to 541 Ma)')
plt.legend(ncol=2,bbox_to_anchor=(0.5, -.02),handles=[Cenozoic_patch,Mesozoic_patch,Paleozoic_late_patch,Paleozoic_early_patch],loc='lower center')
plt.savefig('code_output/0Ma_reconstruction.svg')
plt.show()
t = 250 #reconstruction time
anchor = 1 #anchor plate ID (1 : spin axis)
CEED_blocks_250 = []
pygplates.reconstruct(CEED_blocks, TC17_SHM17_D18_model, CEED_blocks_250, t, anchor)
CEED_land_250 = []
pygplates.reconstruct(CEED_land, TC17_SHM17_D18_model, CEED_land_250, t, anchor)
sutures_active_250 = []
pygplates.reconstruct(sutures_CEED_max_min, TC17_SHM17_D18_model, sutures_active_250, t, anchor)
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
ax.add_patch(patches.Rectangle(xy=[-180, -15], width=360, height=30,
facecolor='#52BE80',
alpha=0.2,
transform=ccrs.PlateCarree()))
plot_reconstructed_feature(ax,CEED_land_250,color='#D7DBDD',linewidth=0.5)
plot_reconstructed_feature(ax,CEED_blocks_250,color='#909497')
plot_reconstructed_feature(ax,sutures_active_250,color='orange',linewidth=2.5)
active_patch = patches.Patch(color='orange',linestyle='-', label='active sutures')
plt.legend(ncol=2,bbox_to_anchor=(0.5, -.02),handles=[active_patch],loc='lower center')
plt.savefig('code_output/250Ma_reconstruction.svg')
plt.show()
t = 300 #reconstruction time
anchor = 1 #anchor plate ID (1 : spin axis)
CEED_blocks_300 = []
pygplates.reconstruct(CEED_blocks, TC17_SHM17_D18_model, CEED_blocks_300, t, anchor)
CEED_land_300 = []
pygplates.reconstruct(CEED_land, TC17_SHM17_D18_model, CEED_land_300, t, anchor)
sutures_active_300 = []
pygplates.reconstruct(sutures_CEED_max_min, TC17_SHM17_D18_model, sutures_active_300, t, anchor)
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
ax.add_patch(patches.Rectangle(xy=[-180, -15], width=360, height=30,
facecolor='#52BE80',
alpha=0.2,
transform=ccrs.PlateCarree()))
plot_reconstructed_feature(ax,CEED_land_300,color='#D7DBDD',linewidth=0.5)
plot_reconstructed_feature(ax,CEED_blocks_300,color='#909497')
plot_reconstructed_feature(ax,sutures_active_300,color='orange',linewidth=2.5)
active_patch = patches.Patch(color='orange',linestyle='-', label='active sutures')
plt.legend(ncol=2,bbox_to_anchor=(0.5, -.02),handles=[active_patch],loc='lower center')
plt.savefig('code_output/300Ma_reconstruction.svg')
plt.show()
t = 445 #reconstruction time
anchor = 1 #anchor plate ID (1 : spin axis)
CEED_blocks_445 = []
pygplates.reconstruct(CEED_blocks, TC17_SHM17_D18_model, CEED_blocks_445, t, anchor)
CEED_land_445 = []
pygplates.reconstruct(CEED_land, TC17_SHM17_D18_model, CEED_land_445, t, anchor)
sutures_active_445 = []
pygplates.reconstruct(sutures_CEED_max_min, TC17_SHM17_D18_model, sutures_active_445, t, anchor)
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
ax.add_patch(patches.Rectangle(xy=[-180, -15], width=360, height=30,
facecolor='#52BE80',
alpha=0.2,
transform=ccrs.PlateCarree()))
plot_reconstructed_feature(ax,CEED_land_445,color='#D7DBDD',linewidth=0.5)
plot_reconstructed_feature(ax,CEED_blocks_445,color='#909497')
plot_reconstructed_feature(ax,sutures_active_445,color='orange',linewidth=2.5)
active_patch = patches.Patch(color='orange',linestyle='-', label='active sutures')
plt.legend(ncol=2,bbox_to_anchor=(0.5, -.02),handles=[active_patch],loc='lower center')
plt.savefig('code_output/445Ma_reconstruction.svg')
plt.show()
These reconstructions are done using the revised Torsvik and Cocks (2017) model (with Swanson-Hysell and Macdonald (2017) and Domeier (2018) adjustments) and the Matthews et al. (2016) model.
time_list_5Myr = range(0,525,5)
anchor = 1
band_width = 10
lat_mins = np.arange(-90 , 90 , band_width)
lat_maxs = np.arange(-90+band_width, 90+band_width, band_width)
lat_mids = lat_mins + (lat_maxs-lat_mins)/2
active_lengths_array_TC17_SHM17_D18_10degree = []
for t in time_list_5Myr:
sutures_active_recon = []
pygplates.reconstruct(sutures_CEED_max_min, TC17_SHM17_D18_model, sutures_active_recon, t, anchor)
active_lengths, active_polylines, active_attributes = get_lengths_in_bands(sutures_active_recon, lat_mins, lat_maxs)
active_lengths_array_TC17_SHM17_D18_10degree.append(active_lengths)
band_width_5degree = 5
lat_mins_5degree = np.arange(-90 , 90 , band_width_5degree)
lat_maxs_5degree = np.arange(-90+band_width_5degree, 90+band_width_5degree, band_width_5degree)
lat_mids_5degree = lat_mins_5degree + (lat_maxs_5degree-lat_mins_5degree)/2.0
active_lengths_array_TC17_SHM17_D18_5degree = []
total_active_length_TC17_SHM17_D18 = []
for t in time_list_5Myr:
sutures_active_recon = []
pygplates.reconstruct(sutures_CEED_max_min, TC17_SHM17_D18_model, sutures_active_recon, t, anchor)
active_lengths, active_polylines, active_attributes = get_lengths_in_bands(sutures_active_recon, lat_mins_5degree, lat_maxs_5degree)
active_lengths_array_TC17_SHM17_D18_5degree.append(active_lengths)
total_active_length_TC17_SHM17_D18.append(np.sum(active_lengths))
time_list_250Ma_5Myr = range(0,255,5)
active_lengths_array_M2016_5degree = []
total_active_length_M2016 = []
for t in time_list_250Ma_5Myr:
sutures_active_recon = []
pygplates.reconstruct(sutures_CEED_max_min, M2016_model, sutures_active_recon, t, anchor)
active_lengths, active_polylines, active_attributes = get_lengths_in_bands(sutures_active_recon, lat_mins_5degree, lat_maxs_5degree)
active_lengths_array_M2016_5degree.append(active_lengths)
total_active_length_M2016.append(np.sum(active_lengths))
length_neg5_5_list_TC17_SHM17_D18 = []
for t in range(0,len(time_list_5Myr)):
length_0_5 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2]
length_neg5_0 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-1]
length_neg5_5 = length_neg5_0 + length_0_5
length_neg5_5_list_TC17_SHM17_D18.append(length_neg5_5)
length_neg10_10_list_TC17_SHM17_D18 = []
for t in range(0,len(time_list_5Myr)):
length_0_5 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2]
length_5_10 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2+1]
length_neg5_0 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-1]
length_neg10_neg5 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-2]
length_neg10_10 = length_neg10_neg5 + length_neg5_0 + length_0_5 + length_5_10
length_neg10_10_list_TC17_SHM17_D18.append(length_neg10_10)
length_neg15_15_list_TC17_SHM17_D18 = []
for t in range(0,len(time_list_5Myr)):
length_0_5 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2]
length_5_10 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2+1]
length_10_15 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2+2]
length_neg5_0 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-1]
length_neg10_neg5 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-2]
length_neg15_neg10 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-3]
length_neg15_15 = length_neg15_neg10 + length_neg10_neg5 + length_neg5_0 + length_0_5 + length_5_10 + length_10_15
length_neg15_15_list_TC17_SHM17_D18.append(length_neg15_15)
length_neg20_20_list_TC17_SHM17_D18 = []
for t in range(0,len(time_list_5Myr)):
length_0_5 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2]
length_5_10 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2+1]
length_10_15 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2+2]
length_15_20 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2+3]
length_neg5_0 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-1]
length_neg10_neg5 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-2]
length_neg15_neg10 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-3]
length_neg20_neg15 = active_lengths_array_TC17_SHM17_D18_5degree[t][len(lat_mids_5degree)/2-4]
length_neg20_20 = length_neg20_neg15 + length_neg15_neg10 + length_neg10_neg5 + length_neg5_0 + length_0_5 + length_5_10 + length_10_15 + length_15_20
length_neg20_20_list_TC17_SHM17_D18.append(length_neg20_20)
length_neg5_5_list_M2016 = []
for t in range(0,len(time_list_250Ma_5Myr)):
length_0_5 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2]
length_neg5_0 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-1]
length_neg5_5 = length_neg5_0 + length_0_5
length_neg5_5_list_M2016.append(length_neg5_5)
length_neg10_10_list_M2016 = []
for t in range(0,len(time_list_250Ma_5Myr)):
length_0_5 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2]
length_5_10 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2+1]
length_neg5_0 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-1]
length_neg10_neg5 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-2]
length_neg10_10 = length_neg10_neg5 + length_neg5_0 + length_0_5 + length_5_10
length_neg10_10_list_M2016.append(length_neg10_10)
length_neg15_15_list_M2016 = []
for t in range(0,len(time_list_250Ma_5Myr)):
length_0_5 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2]
length_5_10 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2+1]
length_10_15 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2+2]
length_neg5_0 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-1]
length_neg10_neg5 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-2]
length_neg15_neg10 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-3]
length_neg15_15 = length_neg15_neg10 + length_neg10_neg5 + length_neg5_0 + length_0_5 + length_5_10 + length_10_15
length_neg15_15_list_M2016.append(length_neg15_15)
length_neg20_20_list_M2016 = []
for t in range(0,len(time_list_250Ma_5Myr)):
length_0_5 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2]
length_5_10 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2+1]
length_10_15 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2+2]
length_15_20 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2+3]
length_neg5_0 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-1]
length_neg10_neg5 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-2]
length_neg15_neg10 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-3]
length_neg20_neg15 = active_lengths_array_M2016_5degree[t][len(lat_mids_5degree)/2-4]
length_neg20_20 = length_neg20_neg15 + length_neg15_neg10 + length_neg10_neg5 + length_neg5_0 + length_0_5 + length_5_10 + length_10_15 + length_15_20
length_neg20_20_list_M2016.append(length_neg20_20)
length_greater40_list_TC17_SHM17_D18 = []
for t in range(0,len(time_list_5Myr)):
length_neg90_neg40 = active_lengths_array_TC17_SHM17_D18_5degree[t][0:10].sum()
length_40_90 = active_lengths_array_TC17_SHM17_D18_5degree[t][-10:].sum()
length_greater40_list_TC17_SHM17_D18.append(length_neg90_neg40 + length_40_90)
plt.figure(figsize=(5,3))
plt.plot(time_list_5Myr,total_active_length_TC17_SHM17_D18,'.-',linewidth=2,label='active suture length')
plt.xlim(520, 0)
plt.xlabel('Age (Ma)')
plt.ylabel('suture length (km)')
plt.title('total suture length')
plt.legend(fontsize=9,loc=2)
plt.show()
plt.figure(figsize=(5,3))
plt.plot(time_list_5Myr,length_neg10_10_list_TC17_SHM17_D18,'.-',
linewidth=2,label='within 10$^\circ$ of equator (TC17_SHM17_D18)')
plt.plot(time_list_250Ma_5Myr,length_neg10_10_list_M2016,'.-',
linewidth=2,label='within 10$^\circ$ of equator (M2016)')
plt.xlim(525, 0)
plt.xlabel('Age (Ma)')
plt.ylabel('suture length (km)')
plt.legend(fontsize=9,loc=2)
plt.show()
plt.figure(figsize=(5,3))
plt.plot(time_list_5Myr,length_neg20_20_list_TC17_SHM17_D18,'.-',
linewidth=2,label='within 20$^\circ$ of equator (TC17_SHM17_D18)')
plt.plot(time_list_250Ma_5Myr,length_neg20_20_list_M2016,'.-',
linewidth=2,label='within 20$^\circ$ of equator (M2016)')
plt.xlim(525, 0)
plt.xlabel('Age (Ma)')
plt.ylabel('suture length (km)')
plt.legend(fontsize=9,loc=2)
plt.show()
plt.figure(figsize=(9,3))
ax1 = plt.subplot(1,2,1)
ax1.plot(time_list_5Myr,total_active_length_TC17_SHM17_D18,'.-',linewidth=2,label='total active suture length',color='k')
plt.xlim(520, 0)
plt.ylim(-4000,36000)
#Hirnantian glaciation (456 to 440 Ma)
Hirnantian = patches.Rectangle((440,-4000),16,4000,edgecolor='none',facecolor='#447CC0')
ax1.add_patch(Hirnantian)
#Devonian glaciation (361-358 Ma)
Devonian = patches.Rectangle((358,-4000),3,4000,edgecolor='none',facecolor='#447CC0')
ax1.add_patch(Devonian)
#Late Paleozoic ice age (337 to 276 Ma)
LPIA = patches.Rectangle((276,-4000),61,4000,edgecolor='none',facecolor='#447CC0')
ax1.add_patch(LPIA)
#Cenozoic glaciation (35 to 0 Ma)
Cenozoic = patches.Rectangle((0,-4000),35,4000,edgecolor='none',facecolor='#447CC0')
ax1.add_patch(Cenozoic)
ax1.hlines(0,0,525)
#ax1.text(535,-1000,'glacial\nintervals',verticalalignment='top',horizontalalignment='right',color='#447CC0')
plt.xlabel('Age (Ma)')
plt.ylabel('suture length (km)')
plt.legend(fontsize=9,loc='upper center')
plt.text(500,30000,'A', color='k',ha ='left',size=20)
ax2 = plt.subplot(1,2,2)
ax2.plot(time_list_5Myr,length_neg10_10_list_TC17_SHM17_D18,'.-',linewidth=2,label='within 10$^\circ$ of equator',color='#481567')
ax2.plot(time_list_5Myr,length_neg15_15_list_TC17_SHM17_D18,'.-',linewidth=2,label='within 15$^\circ$ of equator',color='#33638D')
ax2.plot(time_list_5Myr,length_neg20_20_list_TC17_SHM17_D18,'.-',linewidth=2,label='within 20$^\circ$ of equator',color='#20A387')
plt.xlim(520, 0)
plt.ylim(-2000,18000)
#Hirnantian glaciation (456 to 440 Ma)
Hirnantian = patches.Rectangle((440,-4000),16,4000,edgecolor='none',facecolor='#447CC0')
ax2.add_patch(Hirnantian)
#Devonian glaciation (361-358 Ma)
Devonian = patches.Rectangle((358,-4000),3,4000,edgecolor='none',facecolor='#447CC0')
ax2.add_patch(Devonian)
#Late Paleozoic ice age (337 to 276 Ma)
LPIA = patches.Rectangle((276,-4000),61,4000,edgecolor='none',facecolor='#447CC0')
ax2.add_patch(LPIA)
#Cenozoic glaciation (35 to 0 Ma)
Cenozoic = patches.Rectangle((0,-4000),35,4000,edgecolor='none',facecolor='#447CC0')
ax2.add_patch(Cenozoic)
ax2.hlines(0,0,525)
ax2.text(535,-500,'glacial\nintervals',verticalalignment='top',horizontalalignment='right',color='#447CC0')
plt.xlabel('Age (Ma)')
plt.ylabel('suture length (km)')
plt.legend(fontsize=9,loc='upper center')
plt.text(500,15000,'B', color='k',ha ='left',size=20)
plt.tight_layout()
plt.show()
ice_extent = pd.read_csv('data/ice/IceExtent_Myr.csv')
ice_extent.head()
ice_extent_5Myr = []
ice_extent_5Myr_flip = []
for t in time_list_5Myr:
ice_extent_5Myr.append(ice_extent['LATITUDE FROM POLE'][t])
ice_extent_5Myr_flip.append(90-ice_extent['LATITUDE FROM POLE'][t])
plt.plot(ice_extent['AGE (MA)'],ice_extent['LATITUDE FROM POLE'],label='ice extent (1 Myr sampling)')
plt.plot(time_list_5Myr,ice_extent_5Myr,label='ice extent (5 Myr sampling)')
plt.xlim(525,-1)
plt.legend()
plt.xlabel('Age (Ma)')
plt.ylabel('degrees from pole')
plt.show()
plt.figure(figsize=(5,8))
ax1 = plt.subplot(3,1,1)
ax1.plot(time_list_5Myr,total_active_length_TC17_SHM17_D18,'.-',linewidth=2,label='total suture length',color='k')
plt.xlim(520, -1)
plt.ylim(-1000,30000)
#plt.xlabel('Age (Ma)')
plt.ylabel('active suture length (km)')
plt.legend(fontsize=9,loc='upper center')
plt.text(500,25000,'A', color='k',ha ='left',size=20)
ax1.set_xticklabels([])
ax2 = plt.subplot(3,1,2)
ax2.plot(time_list_5Myr,length_neg10_10_list_TC17_SHM17_D18,'.-',linewidth=2,label='within 10$^\circ$ of equator',color='#481567')
ax2.plot(time_list_5Myr,length_neg15_15_list_TC17_SHM17_D18,'.-',linewidth=2,label='within 15$^\circ$ of equator',color='#33638D')
ax2.plot(time_list_5Myr,length_neg20_20_list_TC17_SHM17_D18,'.-',linewidth=2,label='within 20$^\circ$ of equator',color='#20A387')
plt.xlim(520, -1)
plt.ylim(-500,15000)
plt.ylabel('active suture length (km)')
plt.legend(fontsize=9,loc='upper center')
plt.text(500,12500,'B', color='k',ha ='left',size=20)
ax2.set_xticklabels([])
ax3 = plt.subplot(3,1,3)
#ax3.plot(time_list_5Myr,ice_extent_5Myr_flip)
#ax3.fill_between(time_list_5Myr,ice_extent_5Myr_flip,90)
ax3.plot(ice_extent['AGE (MA)'],90-ice_extent['LATITUDE FROM POLE'])
ax3.fill_between(ice_extent['AGE (MA)'],90-ice_extent['LATITUDE FROM POLE'],90)
plt.xlim(520, -1)
plt.ylim(90.1, 35)
plt.text(500,43,'C', color='k',ha ='left',size=20)
plt.ylabel('ice extent (latitude)')
plt.xlabel('Age (Ma)')
plt.tight_layout()
plt.savefig('code_output/suture_length.pdf')
plt.show()
The tropical LIP area data file that is imported here was generated through the analysis conducted in the LIP_analysis.ipynb notebook within this repository.
LIP_data = pd.read_csv('code_output/tropical_LIP_areas.csv')
LIP_data.sort_values('age (Ma)',inplace=True)
LIP_data.reset_index(drop=True,inplace=True)
LIP_data.head()
plt.figure(figsize=(5,8))
ax1 = plt.subplot(3,1,1)
ax1.plot(LIP_data['age (Ma)'],LIP_data['natural decay total'],label='total (decay scenario)',color='black')
ax1.plot(LIP_data['age (Ma)'],LIP_data['burial total'],label='total (decay + burial scenario)',linestyle='dashed')
plt.ylabel('LIP area (km$^2$)')
plt.legend(fontsize=9,loc='upper right')
plt.text(0.02,0.86,'A', color='k',ha ='left',size=20,transform=ax1.transAxes)
ax1.set_xticklabels([])
plt.xlim(520, -1)
plt.ylim(0, 2.5E7)
ax2 = plt.subplot(3,1,2)
ax2.plot(LIP_data['age (Ma)'],LIP_data['natural decay within 15'],label='within 15$^\circ$ of equator (decay scenario)',color='black')
ax2.plot(LIP_data['age (Ma)'],LIP_data['burial within 15'],linestyle='dashed',label='within 15$^\circ$ of equator (decay + burial scenario)')
plt.xlim(520, -1)
plt.ylim(0, 1.1E7)
plt.ylabel('LIP area (km$^2$)')
plt.legend(fontsize=9,loc='upper right')
plt.text(0.02,0.86,'B', color='k',ha ='left',size=20,transform=ax2.transAxes)
ax2.set_xticklabels([])
ax3 = plt.subplot(3,1,3)
#ax3.plot(time_list_5Myr,ice_extent_5Myr_flip)
#ax3.fill_between(time_list_5Myr,ice_extent_5Myr_flip,90)
ax3.plot(ice_extent['AGE (MA)'],90-ice_extent['LATITUDE FROM POLE'])
ax3.fill_between(ice_extent['AGE (MA)'],90-ice_extent['LATITUDE FROM POLE'],90)
plt.xlim(520, -1)
plt.ylim(90.1, 35)
plt.text(0.02,0.86,'C', color='k',ha ='left',size=20,transform=ax3.transAxes)
plt.ylabel('ice extent (latitude)')
plt.xlabel('Age (Ma)')
plt.tight_layout()
plt.savefig('code_output/LIP_area.pdf')
plt.show()
continental_arc_length = pd.read_csv('data/continental_arcs/Cao2017_arclength.csv')
continental_arc_length.head()
continental_arc_length_5Myr = []
for t in time_list_5Myr:
continental_arc_length_5Myr.append(continental_arc_length['average length (km)'][t])
plt.figure(figsize=(5,3))
ax1 = plt.subplot(1,1,1)
ax1.plot(time_list_5Myr,total_active_length_TC17_SHM17_D18,'.-',linewidth=2,label='active suture length (this study)',color='k')
ax1.plot(time_list_5Myr,continental_arc_length_5Myr,'.-',linewidth=2,label='continental arc length (Cao et al. 2017)',color='#20A387')
plt.xlim(520, -1)
plt.ylim(-1000,45000)
plt.legend(fontsize=9,loc='upper left')
plt.ylabel('length (km)')
plt.xlabel('Age (Ma)')
plt.show()
The dataframe has the following columns:
output_df = pd.DataFrame(columns=['age_Ma'])
output_df['age_Ma'] = time_list_5Myr
output_df['ice_extent'] = ice_extent_5Myr
output_df['total_suture'] = total_active_length_TC17_SHM17_D18
output_df['within_10_suture'] = length_neg10_10_list_TC17_SHM17_D18
output_df['within_15_suture'] = length_neg15_15_list_TC17_SHM17_D18
output_df['within_20_suture'] = length_neg20_20_list_TC17_SHM17_D18
output_df['greater_10_suture'] = output_df['total_suture'] - output_df['within_10_suture']
output_df['greater_15_suture'] = output_df['total_suture'] - output_df['within_15_suture']
output_df['greater_20_suture'] = output_df['total_suture'] - output_df['within_20_suture']
output_df['greater_40_suture'] = length_greater40_list_TC17_SHM17_D18
output_df['total_LIP_decay'] = LIP_data['natural decay total']
output_df['total_LIP_decay_burial'] = LIP_data['burial total']
output_df['within_15_LIP_decay'] = LIP_data['natural decay within 15']
output_df['within_15_LIP_decay_burial'] = LIP_data['burial within 15']
output_df['outside_15_LIP_decay'] = LIP_data['natural decay outside 15']
output_df['outside_15_LIP_decay_burial'] = LIP_data['burial outside 15']
output_df['continental_arc_Cao2017'] = continental_arc_length_5Myr
output_df.to_csv('code_output/ice_LIP_suture_lengths.csv',index=False)
output_df.head()
plt.figure(figsize=(5,12))
ax1 = plt.subplot(3,1,1)
#plot within_20 + greater_20
ax1.plot(time_list_5Myr,output_df['greater_20_suture']+output_df['within_20_suture'],color='k',label='total suture length')
#plot fill between total suture and within_20 suture length
ax1.fill_between(time_list_5Myr,output_df['within_15_suture'],output_df['total_suture'],label='>15 suture length')
#plot fill between within_20 suture length and 0
ax1.fill_between(time_list_5Myr,output_df['within_15_suture'],0,label='<15 suture length')
plt.xlim(520, -1)
plt.ylabel('length (km)')
plt.xlabel('Age (Ma)')
plt.legend()
plt.show()
plt.scatter(output_df['greater_40_suture'],output_df['ice_extent'],label='>40 degrees suture length')
plt.scatter(output_df['within_15_suture'],output_df['ice_extent'],label='<15 degrees suture length')
plt.xlabel('suture length (km)')
plt.ylabel('degrees from pole')
plt.xlim(-500,18000)
plt.ylim(-5,65)
plt.legend(loc='lower right')
plt.show()
from scipy.stats.stats import pearsonr
print('pearsonr between <15 sutures length and ice extent')
print(pearsonr(output_df['within_15_suture'],ice_extent_5Myr))
print('')
print('pearsonr between >15 sutures length and ice extent')
print(pearsonr(output_df['greater_15_suture'],ice_extent_5Myr))
print('')
print('pearsonr between <20 sutures length and ice extent')
print(pearsonr(output_df['within_20_suture'],ice_extent_5Myr))
print('')
print('pearsonr between >20 sutures length and ice extent')
print(pearsonr(output_df['greater_20_suture'],ice_extent_5Myr))
print('')
print('pearsonr between >40 sutures length and ice extent')
print(pearsonr(output_df['greater_40_suture'],ice_extent_5Myr))
print('')
print('pearsonr between total sutures length and ice extent')
print(pearsonr(output_df['total_suture'],ice_extent_5Myr))
print('')
print('pearsonr between total LIP area (decay) and ice extent')
print(pearsonr(output_df['total_LIP_decay'],ice_extent_5Myr))
print('')
print('pearsonr between <15 LIP area (decay) and ice extent')
print(pearsonr(output_df['within_15_LIP_decay'],ice_extent_5Myr))
print('')
print('pearsonr between total LIP area (burial + decay) and ice extent')
print(pearsonr(output_df['total_LIP_decay_burial'],ice_extent_5Myr))
print('')
print('pearsonr between <15 LIP area (burial + decay) and ice extent')
print(pearsonr(output_df['within_15_LIP_decay_burial'],ice_extent_5Myr))
print('')
print('pearsonr between continental arc length and ice extent')
print(pearsonr(continental_arc_length_5Myr,ice_extent_5Myr))
fig, axes = plt.subplots(4,3, sharex=True, sharey=True,figsize=(9,9))
plt.xlim(-90,90)
plt.xticks(np.arange(-90, 120, step=30))
for i, ax in enumerate(axes.flatten()):
if i< 11:
ax.bar(lat_mids, active_lengths_array_TC17_SHM17_D18_10degree[i*10],width=10,color='grey',align='center')
tropics = patches.Rectangle((-15,0),30,9000,edgecolor='none',alpha=0.2,facecolor='#52BE80')
ax.add_patch(tropics)
plt.ylim(0,9000)
if float(time_list_5Myr[i*10]) < 35:
ax.text(-80,8000,str(time_list_5Myr[i*10])+' Ma', color='blue')
elif float(time_list_5Myr[i*10]) < 276:
ax.text(-80,8000,str(time_list_5Myr[i*10])+' Ma', color='red')
elif float(time_list_5Myr[i*10]) < 337:
ax.text(-80,8000,str(time_list_5Myr[i*10])+' Ma', color='blue')
elif float(time_list_5Myr[i*10]) < 440:
ax.text(-80,8000,str(time_list_5Myr[i*10])+' Ma', color='red')
elif float(time_list_5Myr[i*10]) < 456:
ax.text(-80,8000,str(time_list_5Myr[i*10])+' Ma', color='blue')
elif float(time_list_5Myr[i*10]) < 540:
ax.text(-80,8000,str(time_list_5Myr[i*10])+' Ma', color='red')
axes.flatten()[0].set_ylabel('suture length (km)')
axes.flatten()[3].set_ylabel('suture length (km)')
axes.flatten()[6].set_ylabel('suture length (km)')
axes.flatten()[9].set_ylabel('suture length (km)')
#axes.flatten()[12].set_ylabel('suture length (km)')
axes.flatten()[-3].set_xlabel('latitude')
axes.flatten()[-2].set_xlabel('latitude')
axes.flatten()[-1].set_xlabel('latitude')
title = fig.suptitle("reconstructed active ophiolite-bearing sutures (TC2017 with S-HM2017 and D2018 modifications)", fontsize=14, horizontalalignment='center')
title.set_y(1)
title.set_x(0.525)
plt.tight_layout()
plt.savefig('code_output/active_suture_bar_plot.pdf')
plt.show()
This code creates snapshots of reconstructions in 5 million year steps over the past 500 million years. These images are combined in the suture_GIF_creation.ipynb notebook to make animations.
anchor = 1 #anchor plate ID (1 : spin axis)
recon_image_list = []
for t in range(0,525,5):
CEED_blocks_recon = []
pygplates.reconstruct(CEED_blocks, TC17_SHM17_D18_model, CEED_blocks_recon, t, anchor)
sutures_active_recon = []
pygplates.reconstruct(sutures_CEED_max_min, TC17_SHM17_D18_model, sutures_active_recon, t, anchor)
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
ax.add_patch(patches.Rectangle(xy=[-180, -15], width=360, height=30,
facecolor='#52BE80',
alpha=0.2,
transform=ccrs.PlateCarree()))
plot_reconstructed_feature(ax,CEED_blocks_recon)
plot_reconstructed_feature(ax,sutures_active_recon,color='orange',linewidth=2.5)
plt.title(str(t) + ' Ma reconstruction\n (active ophiolite-bearing sutures are orange)')
image_name = 'code_output/paleogeo_snapshots/recon' + str(t) + '.png'
recon_image_list.append(image_name)
plt.savefig(image_name)
plt.close()
What is shown in these images is a schematic extent of ice that extends as far as the maximum extent in the compilation. These images are not a representation of ice sheets through time and are developed to be broadly illustrative.
output_df_ageindex = output_df.set_index('age_Ma')
for t in range(0,525,5):
CEED_blocks_recon = []
pygplates.reconstruct(CEED_blocks, TC17_SHM17_D18_model, CEED_blocks_recon, t, anchor)
sutures_active_recon = []
pygplates.reconstruct(sutures_CEED_max_min, TC17_SHM17_D18_model, sutures_active_recon, t, anchor)
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection=ccrs.Robinson()))
ax.set_global()
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
ax.add_patch(patches.Rectangle(xy=[-180, -15], width=360, height=30,
facecolor='#52BE80',
alpha=0.2,
transform=ccrs.PlateCarree()))
if t<10:
ax.add_patch(patches.Rectangle(xy=[-180, 90-output_df_ageindex.ice_extent[t]], width=360, height=output_df_ageindex.ice_extent[t],
facecolor='#37A4DA',
alpha=0.4,
transform=ccrs.PlateCarree()))
ax.add_patch(patches.Rectangle(xy=[-180, -90], width=360, height=output_df_ageindex.ice_extent[t],
facecolor='#37A4DA',
alpha=0.4,
transform=ccrs.PlateCarree()))
if t>=10:
ax.add_patch(patches.Rectangle(xy=[-180, -90], width=360, height=output_df_ageindex.ice_extent[t],
facecolor='#37A4DA',
alpha=0.4,
transform=ccrs.PlateCarree()))
plot_reconstructed_feature(ax,CEED_blocks_recon)
plot_reconstructed_feature(ax,sutures_active_recon,color='orange',linewidth=2.5)
plt.title(str(t) + ' Ma reconstruction\n (active ophiolite-bearing sutures are orange)\n (schematic maximum ice extent in blue)')
image_name = 'code_output/paleogeo_snapshots/recon_ice' + str(t) + '.png'
recon_image_list.append(image_name)
plt.savefig(image_name)
plt.close()
For the sake of documentation, the functions imported at the start of the notebook from the recon_tools.py functions that were developed for these analyses are shown below.
import math
from numpy.core.umath_tests import inner1d
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import numpy as np
from shapely.geometry.polygon import Polygon
import pygplates
import pandas as pd
########## BASIC CALCULATIONS ##########
def lat_lon_2_cart(lat, lon):
"""
Convert lat/lon coordinates to cartesian.
Parameters
----------
lat : array-like
latitude (-90 to 90)
lon : array-like
longitude (-180 to 180)
Returns
-------
cart : tuple
(x, y, z) cartesian coordinates on the unit sphere
"""
# convert to radians
lat = math.radians(lat)
lon = math.radians(lon)
# calculations
x = math.cos(lon) * math.cos(lat)
y = math.sin(lon) * math.cos(lat)
z = math.sin(lat)
cart = (x, y, z)
return cart
def cart_2_lat_lon(cart):
"""
Convert cartesian coordinates to lat/lon.
Parameters
----------
cart : tuple
(x, y, z) cartesian coordinates on the unit sphere
Returns
-------
lat : array-like
latitude (-90 to 90)
lon : array-like
longitude (-180 to 180)
"""
# calculations
lon = math.atan2(cart[1], cart[0])
lat = math.atan2(cart[2], math.sqrt(cart[0]**2 + cart[1]**2))
# convert to degrees
lat = math.degrees(lat)
lon = math.degrees(lon)
return lat, lon
def fast_cross(a, b):
"""
3D matrix cross multiplication.
source: http://ssb.stsci.edu/doc/stsci_python_x/stsci.sphere.doc/html/_modules/stsci/sphere/great_circle_arc.html
Parameters
----------
a, b : numpy arrays
matrices to be cross multiplied
Returns
-------
cp : numpy array
cross product
"""
cp = np.empty(np.broadcast(a, b).shape)
aT = a.T
bT = b.T
cpT = cp.T
cpT[0] = aT[1]*bT[2] - aT[2]*bT[1]
cpT[1] = aT[2]*bT[0] - aT[0]*bT[2]
cpT[2] = aT[0]*bT[1] - aT[1]*bT[0]
return cp
def cross_and_normalize(A, B):
"""
3D matrix cross multiplication and normalized.
source: http://ssb.stsci.edu/doc/stsci_python_x/stsci.sphere.doc/html/_modules/stsci/sphere/great_circle_arc.html
Parameters
----------
A, B : numpy arrays
matrices to be cross multiplied
Returns
-------
TN : numpy array
normalized cross product
"""
T = fast_cross(A, B)
# normalization
l = np.sqrt(np.sum(T ** 2, axis=-1))
l = np.expand_dims(l, 2)
# might get some divide-by-zeros, but we don't care
with np.errstate(invalid='ignore'):
TN = T / l
return TN
def intersection(A, B, C, D):
"""
Point of intersection between two great circle arcs.
source: http://ssb.stsci.edu/doc/stsci_python_x/stsci.sphere.doc/html/_modules/stsci/sphere/great_circle_arc.html
Parameters
----------
A, B : (*x*, *y*, *z*) triples or Nx3 arrays of triples
Endpoints of the first great circle arc.
C, D : (*x*, *y*, *z*) triples or Nx3 arrays of triples
Endpoints of the second great circle arc.
Returns
-------
T : (*x*, *y*, *z*) triples or Nx3 arrays of triples
If the given arcs intersect, the intersection is returned. If the arcs do not intersect,
the triple is set to all NaNs.
"""
A = np.asanyarray(A)
B = np.asanyarray(B)
C = np.asanyarray(C)
D = np.asanyarray(D)
A, B = np.broadcast_arrays(A, B)
C, D = np.broadcast_arrays(C, D)
ABX = fast_cross(A, B)
CDX = fast_cross(C, D)
T = cross_and_normalize(ABX, CDX)
T_ndim = len(T.shape)
if T_ndim > 1:
s = np.zeros(T.shape[0])
else:
s = np.zeros(1)
s += np.sign(inner1d(fast_cross(ABX, A), T))
s += np.sign(inner1d(fast_cross(B, ABX), T))
s += np.sign(inner1d(fast_cross(CDX, C), T))
s += np.sign(inner1d(fast_cross(D, CDX), T))
if T_ndim > 1:
s = np.expand_dims(s, 2)
cross = np.where(s == -4, -T, np.where(s == 4, T, np.nan))
# If they share a common point, it's not an intersection. This
# gets around some rounding-error/numerical problems with the
# above.
equals = (np.all(A == C, axis=-1) |
np.all(A == D, axis=-1) |
np.all(B == C, axis=-1) |
np.all(B == D, axis=-1))
equals = np.expand_dims(equals, 2)
T = np.where(equals, np.nan, cross)
return T
########## ZONAL CALCULATIONS ##########
def check_polygon_in_band(polygon, lat_min, lat_max):
"""
Check to see whether any part of a given polygon is inside a given latitude band.
Parameters
----------
polygon : pygpplates polygon
lat_min : float
the minimum latitude of the latitude band
lat_max : float
the maximum latitude of the latitude band
Returns
-------
in_band : boolean
True if inside, False if outside
"""
# pull out lat/lon vertices
lat_lon_array = polygon.to_lat_lon_array()
lats = lat_lon_array[:,0]
# check to see if any part of the polygon falls into our latitude band
in_band = False
for j in range(len(lats)):
if lats[j]>lat_min and lats[j]<lat_max:
in_band = True
break
return in_band
def get_area_in_band(polygon, lat_min, lat_max):
"""
Calculate the area of a given polygon inside a given latitude band.
Parameters
----------
polygon : pygpplates polygon
lat_min : float
the minimum latitude of the latitude band
lat_max : float
the maximum latitude of the latitude band
Returns
-------
area : float
the area of the polygon within the latitude band (in km^2)
band_polygon : pygplates polygon
with the parts outside of the latitude band removed
"""
# pull out lat/lon vertices
lat_lon_array = polygon.to_lat_lon_array()
lats = lat_lon_array[:,0]
lons = lat_lon_array[:,1]
# storage lists
bookmarks = []
lat_add_list = []
lon_add_list = []
# iterate through the points
for i in range(1,len(lats)):
top_cross = False
bot_cross = False
# case where we move into the band from above
if lats[i-1]>lat_max and lats[i]<lat_max:
top_cross = True
# case where we move out of the band from below
if lats[i-1]<lat_max and lats[i]>lat_max:
top_cross = True
# case where we move out of the band from above
if lats[i-1]>lat_min and lats[i]<lat_min:
bot_cross = True
# case where we move into the band from below
if lats[i-1]<lat_min and lats[i]>lat_min:
bot_cross = True
# do calculations if we cross
if top_cross or bot_cross:
# convert the endpoints of the polygon segment into cartesian
A = lat_lon_2_cart(lats[i-1], lons[i-1])
B = lat_lon_2_cart(lats[i] , lons[i])
# get the intersection point (for the top and bottom cases), and convert back to lat/lon
if top_cross:
C_top = lat_lon_2_cart(lat_max, min([lons[i],lons[i-1]]))
D_top = lat_lon_2_cart(lat_max, max([lons[i],lons[i-1]]))
T = intersection(A, B, C_top, D_top)
else:
C_bot = lat_lon_2_cart(lat_min, min([lons[i],lons[i-1]]))
D_bot = lat_lon_2_cart(lat_min, max([lons[i],lons[i-1]]))
T = intersection(A, B, C_bot, D_bot)
lat_add, lon_add = cart_2_lat_lon(T)
# add to the storage lists
bookmarks.append(i)
lat_add_list.append(lat_add)
lon_add_list.append(lon_add)
# now insert the stored values into the original arrays
new_lats = np.insert(lats, bookmarks, lat_add_list)
new_lons = np.insert(lons, bookmarks, lon_add_list)
# only keep values below the maximum latitude (with small buffer)
mask = np.less(new_lats, lat_max+0.1)
new_lats = new_lats[mask]
new_lons = new_lons[mask]
# only keep values above the minimum latitude
mask = np.greater(new_lats, lat_min-0.1)
new_lats = new_lats[mask]
new_lons = new_lons[mask]
# create a Polygon, if we are left with enough points
if len(new_lats) >= 3:
band_polygon = pygplates.PolygonOnSphere(zip(new_lats,new_lons))
# get the area in km2
area = band_polygon.get_area() * 6371.009**2
# if we don't...
else:
area = 0
band_polygon = None
return area, band_polygon
def get_areas_in_bands(reconstructed_feature_geometries, lat_mins, lat_maxs):
"""
Get the area of all features in each latitude band.
Parameters
----------
reconstructed_feature_geometries : list
list of reconstructed features
lat_mins : array-like
array-like of latitude minimums
lat_maxs : array_like
array_like of latitude maximums
Returns
-------
areas : array
list of total area in each latitude band
area_polygons : list
list of all polygons for which areas were calculated
"""
# storage vectors
areas = np.array([])
area_polygons = []
# iterate over each latitude band
for i in range(len(lat_mins)):
accumulated_area = 0
# iterate over each polygon
for j in range(len(reconstructed_feature_geometries)):
current_polygon = reconstructed_feature_geometries[j].get_reconstructed_geometry()
# check if the polygon is in the band
in_band = check_polygon_in_band(current_polygon, lat_mins[i], lat_maxs[i])
if in_band:
# do the calculation
area, band_polygon = get_area_in_band(current_polygon, lat_mins[i], lat_maxs[i])
# store results
accumulated_area = accumulated_area + area
area_polygons.append(band_polygon)
# store total area for the band
areas = np.append(areas, accumulated_area)
return areas, area_polygons
def get_length_in_band(polyline, lat_min, lat_max):
"""
Calculate the length of a given polyline inside a given latitude band.
Parameters
----------
polyline : pygpplates polyline
lat_min : float
the minimum latitude of the latitude band
lat_max : float
the maximum latitude of the latitude band
Returns
-------
length : float
the length of the polyline within the latitude band (in km)
band_polylines : list
pygplates polylines with the parts outside of the latitude band removed
"""
# pull out lat/lon vertices
lat_lon_array = polyline.to_lat_lon_array()
lats = lat_lon_array[:,0]
lons = lat_lon_array[:,1]
# storage lists
bookmarks = []
lat_add_list = []
lon_add_list = []
# iterate through the points
for i in range(1,len(lats)):
top_cross = False
bot_cross = False
# case where we move into the band from above
if lats[i-1]>lat_max and lats[i]<lat_max:
top_cross = True
# case where we move out of the band from below
if lats[i-1]<lat_max and lats[i]>lat_max:
top_cross = True
# case where we move out of the band from above
if lats[i-1]>lat_min and lats[i]<lat_min:
bot_cross = True
# case where we move into the band from below
if lats[i-1]<lat_min and lats[i]>lat_min:
bot_cross = True
# do calculations if we cross
if top_cross or bot_cross:
# convert the endpoints of the polygon segment into cartesian
A = lat_lon_2_cart(lats[i-1], lons[i-1])
B = lat_lon_2_cart(lats[i] , lons[i])
# get the intersection point (for the top and bottom cases), and convert back to lat/lon
if top_cross:
C_top = lat_lon_2_cart(lat_max, min([lons[i],lons[i-1]]))
D_top = lat_lon_2_cart(lat_max, max([lons[i],lons[i-1]]))
T = intersection(A, B, C_top, D_top)
else:
C_bot = lat_lon_2_cart(lat_min, min([lons[i],lons[i-1]]))
D_bot = lat_lon_2_cart(lat_min, max([lons[i],lons[i-1]]))
T = intersection(A, B, C_bot, D_bot)
lat_add, lon_add = cart_2_lat_lon(T)
# add to the storage lists
bookmarks.append(i)
lat_add_list.append(lat_add)
lon_add_list.append(lon_add)
# now insert the stored values into the original arrays
new_lats = np.insert(lats, bookmarks, lat_add_list)
new_lons = np.insert(lons, bookmarks, lon_add_list)
# mask points above and below the latitude band (with small buffer)
top_mask = np.less(new_lats, lat_max+0.1)
bot_mask = np.greater(new_lats, lat_min-0.1)
mask = top_mask & bot_mask
# initiate final output variabiles
length = 0
band_polylines = []
# initiate lat lon arrays for the current polyline that is being split out
split_lats = np.array([])
split_lons = np.array([])
hold = False
# iterate through the mask
for i in range(len(mask)):
if mask[i]:
# begin holding, and store points
hold = True
split_lats = np.append(split_lats, new_lats[i])
split_lons = np.append(split_lons, new_lons[i])
# i.e. the first time we reach a point outside of the band
elif hold:
# create a polyline, if we are left with enough points
if len(split_lats) >= 2:
band_polyline = pygplates.PolylineOnSphere(zip(split_lats,split_lons))
# get the length in km
length_split = band_polyline.get_arc_length() * pygplates.Earth.mean_radius_in_kms
# store
length = length + length_split
band_polylines.append(band_polyline)
# reset stuff
split_lats = np.array([])
split_lons = np.array([])
hold = False
# if the last point was within the band
if hold:
# create a polyline, if we are left with enough points
if len(split_lats) >= 2:
band_polyline = pygplates.PolylineOnSphere(zip(split_lats,split_lons))
# get the length in km
length_split = band_polyline.get_arc_length() * pygplates.Earth.mean_radius_in_kms
# store
length = length + length_split
band_polylines.append(band_polyline)
return length, band_polylines
def get_lengths_in_bands(reconstructed_feature_geometries, lat_mins, lat_maxs):
"""
Get the area of all features in each latitude band.
Parameters
----------
reconstructed_feature_geometries : list
list of reconstructed features
lat_mins : array-like
array-like of latitude minimums
lat_maxs : array_like
array_like of latitude maximums
Returns
-------
accumulated_lengths : array
array of total area in each latitude band
length_polylines : polylines
list of all polylines for which areas were calculated
polyline_attributes : Pandas Dataframe
Dataframe with names, lengths, latitudes and recon_times for each polyline
"""
# storage vectors
lengths = np.array([])
accumulated_lengths = np.array([])
length_polylines = []
names = []
polyline_lat_min = []
polyline_lat_max = []
recon_times = []
# iterate over each latitude band
for i in range(len(lat_mins)):
accumulated_length = 0
# iterate over each polygon
for j in range(len(reconstructed_feature_geometries)):
current_polyline = reconstructed_feature_geometries[j].get_reconstructed_geometry()
# check if the polygon is in the band
in_band = check_polygon_in_band(current_polyline, lat_mins[i], lat_maxs[i])
if in_band:
# do the calculation
length, band_polyline = get_length_in_band(current_polyline, lat_mins[i], lat_maxs[i])
name = reconstructed_feature_geometries[j].get_feature().get_name()
recon_time = reconstructed_feature_geometries[j].get_reconstruction_time()
# store results
accumulated_length = accumulated_length + length
length_polylines.append(band_polyline)
names.append(name)
lengths = np.append(lengths, length)
polyline_lat_min.append(lat_mins[i])
polyline_lat_max.append(lat_maxs[i])
recon_times.append(recon_time)
# store total area for the band
accumulated_lengths = np.append(accumulated_lengths, accumulated_length)
polyline_attributes = pd.DataFrame({'name':names,'length':lengths,'polyline_lat_min':polyline_lat_min,'polyline_lat_max':polyline_lat_max,'recon_time':recon_times})
return accumulated_lengths, length_polylines, polyline_attributes
########## PLOTTING FUNCTIONS ##########
def plot_feature(ax,feature,color='grey',linewidth=1):
for n in range(len(feature)):
# pull out lat/lon vertices
lat_lon_array = feature[n].to_lat_lon_array()
lats = lat_lon_array[:,0]
lons = lat_lon_array[:,1]
ax.plot(lons,lats, transform=ccrs.Geodetic(), color=color, linewidth=linewidth)
def plot_reconstructed_feature(ax,reconstructed_feature,color='grey',linewidth=1):
for n in range(len(reconstructed_feature)):
# pull out lat/lon vertices
lat_lon_array = reconstructed_feature[n].get_reconstructed_geometry().to_lat_lon_array()
lats = lat_lon_array[:,0]
lons = lat_lon_array[:,1]
ax.plot(lons,lats, transform=ccrs.Geodetic(), color=color, linewidth=linewidth)
def plot_reconstruction(reconstructed_feature_geometries_list, color_list, lon_0):
"""
Plot a global reconstruction from pygplates.
Parameters
----------
reconstructed_feature_geometries_list : list of lists
list of lists of reconstructed features
color_list : list of colors
list of matplotlib color for geometries
lon_0 : float
the central longitude for viewing
Returns
-------
None.
"""
# initialize map
fig = plt.figure(figsize=(12,10))
ax = plt.subplot(1, 1, 1, projection=ccrs.Robinson(central_longitude=lon_0))
ax.set_title(str(reconstructed_feature_geometries_list[0][0].get_reconstruction_time()) + ' Ma')
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
# loop over each reconstructed geometry list
for i in range(len(reconstructed_feature_geometries_list)):
# loop over each reconstructed geometry
for j in range(len(reconstructed_feature_geometries_list[i])):
# pull out lat/lon vertices
lat_lon_array = reconstructed_feature_geometries_list[i][j].get_reconstructed_geometry().to_lat_lon_array()
lats = lat_lon_array[:,0]
lons = lat_lon_array[:,1]
# zip the result
poly = Polygon(zip(lons, lats))
# add the polygon to the map
ax.add_geometries([poly], ccrs.PlateCarree(), facecolor=color_list[i], edgecolor='k', alpha=0.5)
#ax.plot(lons,lats,transform=ccrs.Geodetic(),color=color_list[i],linewidth=1)
plt.show()
def plot_polygons(polygon_list, color, lon_0):
"""
Plot pygplates polygons.
Parameters
----------
polygon_list : list
list of pygplates polygons
color : color
matplotlib color for geometries
lon_0 : float
the central longitude for viewing
Returns
-------
None.
"""
# initialize map
fig = plt.figure(figsize=(12,10))
ax = plt.subplot(1, 1, 1, projection=ccrs.Robinson(central_longitude=lon_0))
ax.gridlines(xlocs=np.arange(-180,181,60),ylocs=np.arange(-90,91,30),linestyle='--')
# loop over each polygon
for i in range(len(polygon_list)):
if polygon_list[i] != None:
# pull out lat/lon vertices
lat_lon_array = polygon_list[i].to_lat_lon_array()
lats = lat_lon_array[:,0]
lons = lat_lon_array[:,1]
# zip the result
poly = Polygon(zip(lons, lats))
# add the polygon to the map
ax.add_geometries([poly], ccrs.PlateCarree(), facecolor=color, edgecolor='k', alpha=0.5)
#ax.plot(lons,lats,transform=ccrs.Geodetic(),color=color,linewidth=1)
plt.show()